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Abstract 

We examine several numerical techniques for the calculation of the dynamics of quantum systems. In particular, we 
single out an iterative method which is based on expanding the time evolution operator into a finite series of Chebyshev 
polynomials. The Chebyshev approach benefits from two advantages over the standard time-integration Crank- Nicholson 
scheme: speedup and efficiency. Potential competitors are semiclassical methods such as the Wigner-Moyal or quantum 
tomographic approaches. We outline the basic concepts of these techniques and benchmark their performance against 
the Chebyshev approach by monitoring the time evolution of a Gaussian wave packet in restricted one-dimensional (ID) 
geometries. Thereby the focus is on tunnelling processes and the motion in anharmonic potentials. Finally we apply the 
prominent Chebyshev technique to two highly non-trivial problems of current interest: (i) the injection of a particle in 
a disordered 2D graphene nanoribbon and (ii) the spatiotemporal evolution of polaron states in finite quantum systems. 
Here, depending on the disorder/electron-phonon coupling strength and the device dimensions, we observe transmission 
or localisation of the matter wave. 
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1. Introduction 

Quantum statistical physics, such as condensed matter 
or plasma physics, but also quantum chemistry, heavily 
depends on effective numerical methods for solving com- 
plex few-particle and many particle problems. Implement- 
ing suitable theoretical concepts for their description on 
modern supercomputer architectures, nowadays computa- 
tional physics constitutes, besides experiment and theory, 
the third pillar of contemporary physics Q- Numerical 
techniques become especially important for strongly cor- 
related systems where analytical approaches largely fail. 
This is due to the absence of small (coupling) parame- 
ters or, more general, because the relevant energy scales 
are not well separated, both preventing the application of 
standard perturbative schemes. 

Common to any numerical method in quantum physics 
is the requirement to represent the states and operators 
describing the physical system in the Hilbert space in a 
form that is suited for computations. Then, working with 
a discrete basis of the Hilbert space, the computational 
challenge is the solution of an eigenvalue problem for huge 
(sparse) matrices. For most physical systems the dimen- 
sion of the Hilbert space is much too large in order to per- 
form a full exact diagonalisation of the related Hamilton 
matrix. Fortunately some quantities of interest depend on 
the properties of the ground state or a few excited states 
only, and therefore may be studied by iterative Krylov 
space techniques such as Lanczos diagonalisation The 
quantum dynamics or long time behaviour of correlated 
systems, however, require, in principle, the knowledge of 



all eigenstates. 

The theoretical investigation of quantum dynamics was 
triggered in recent years by the vast progress of the ex- 
perimental techniques. Nowadays femtosecond laser spec- 
troscopy, for instance, allows for a precise analysis of quan- 
tum dynamical processes with extreme time resolution. 
Direct time integration of the Schrodinger equation at the 
cost of a full diagonalisation of the system's Hamiltonian 
(including the coupling to external fields) is impractical in 
such cases because of its computational complexity. 

The aim of this work is to propose a very efficient 
Chebyshev-based algorithm that allows calculating the dy- 
namics of a quantum system numerically exactly, also for 
relatively long times, and therefore overcomes the above 
mentioned problem at least partially. In order to demon- 
strate the power of our iterative Chebyshev expansion ap- 
proach, we compare the accuracy and computational costs 
of certain model calculations with those emerging by the 
use of the more the standard Crank-Nicholson, Wigner- 
Moyal and quantum tomography methods. We start by 
presenting the basic ideas of the iterative (Chebyshev, Crank- 
Nicholson; Sect. 12. lj ) and semiclassical (Wigner-Moyal, to- 
mographic; Sect. 12.3ft techniques. Afterwards, in Sect. 
we consider three different problems of increasing com- 
plexity: (i) the motion of a Gaussian wave packet in a ID 
geometry (Sec. l3.lj >. (ii) the evolution of a particle in a dis- 
ordered 2D graphene nanoribbon (Sec. 13. 2| ). and (iii) the 
spatiotemporal evolution of polaron states in finite quan- 
tum systems (Sec. 13.31) . Our conclusions will be presented 
in Sect.H 
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2. Time evolution of quantum systems 

The time evolution of a quantum state \ip) is described 
by the Schrodinger equation 

m^- t m))=H\m)' (i) 

If the Hamilton operator H does not explicitly depend 
on time t we can formally integrate this equation and ex- 
press the dynamics of a given quantum state \ip(to)) in 
terms of the time evolution operator U(t,to) as \ip(t)) = 
U(t,t )\ip(t )), where U(t,t ) = e -iH(*-*o)/S, Exploiting 
that U(t, tg) is diagonal in the eigenbasis of the Hamilto- 
nian, we can directly determine the dynamics of the quan- 
tum system. 

2.1. Direct method 

For systems with moderate Hilbert space dimensions, a 
full diagonalisation of the Hamiltonian permits expression 
of the quantum dynamics of an initial state \ip(to)) &s 
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Here \n) are the (time independent) eigenstates of the sys- 
tem and E n the corresponding eigenenergies. In this way 
the decomposition of an initial state into a linear combi- 
nation of eigenstates allows for an exact calculation of the 
quantum state at arbitrary times. As soon as the physical 
description of the system requires a larger Hilbert space 
dimension, however, this direct calculation is no longer 
feasible and we have to resort to alternative approaches. 

2.2. Iterative methods 

Crank-Nicholson scheme. One of the standard meth- 
ods to calculate the quantum time evolution iteratively is 
the Crank- Nicholson algorithm 18]. Dividing [to,t] into 
S subintervals 5t, the quantum state evolves for each it- 
erative time step St from t s to t s+ i = t s + 5t according 
to 
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There are two limitations to this scheme. First, in addition 
to the matrix vector multiplication (MVM) on the right 
hand side of ((3]), each iteration requires the solution of a 
linear system of equations to obtain \tp(t s +i)). Despite the 
availability of many powerful methods for the solution of 
large (sparse) linear equation systems, this task remains 
the most time consuming part of the algorithm. Using 
iterative methods for the solution of the linear equation 
system, the attainable Hilbert space dimensions increase 
substantially as compared to direct methods. Second, the 
Crank-Nicholson algorithm is accurate only to order (St) 2 , 
which severely restricts the maximum usable iterative time 
step. 



Chebyshev scheme. Both limitations can be overcome 
by an approach where we expand the time evolution op- 
erator U(t,to) = U(t — to) = U(At) into a finite series 
of first-kind Chebyshev polynomials of order k: Tk(x) = 
cos(fcarccos(a;)). We then obtain 0,0, ED] 
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(4) 

Prior to the expansion the Hamiltonian has to be shifted 
and rescaled such that the spectrum of H = (H — b)/a is 
within the definition interval of the Chebyshev polynomi- 
als, [—1,1]. The parameters a and b are calculated from 
the extreme eigenvalues of H as b — |(£ max + E m i n ) and 
a = l(E max -E min + e). Here we introduced e = ot(E max - 
-Emin) to ensure the rescaled spectrum \E\ < 1/(1 + a) lies 
well inside [—1,1]. In practice, we use a = 0.01. Note that 
the Chebyshev expansion also applies to systems with un- 
bounded spectra. In those cases we truncate the infinite 
Hilbert space to a finite dimension by restricting the model 
on a discrete space grid or using an energy cutoff. In this 
way we ensure the finiteness of the extreme eigenvalues. 
In the expansion coefficients c k are given by 
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(5) 

( Jfc denotes the fc-th order Bessel function of the first kind). 

To calculate the evolution of a state \ip(to)) horn one 
time grid point to the adjacent one, \ip(t)) = U(At)\ip(to)), 
we have to accumulate the (^--weighted vectors \v k ) = 
T k (H)\ip(to)). Since the coefficients c k (aAt/%) depend on 
the time step but not on time explicitly, we need to cal- 
culate them only once. The vectors \v k ) can be computed 
iteratively exploiting the recurrence relation of the Cheby- 
shev polynomials, 



\v k+1 ) = 2H\v k ) - \Vk-i) 



(6) 



with \v\) — H\vq) and |uo) = \ip(ta)). Evolving the wave 
function from one time step to the next requires M MVMs 
of a given complex vector with the (sparse) Hamilton ma- 
trix of dimension N and the summation of the resulting 
vectors after an appropriate rescaling. The Chebyshev ex- 
pansion may also be applied to systems with time depen- 
dent Hamiltonians, but there the time variation H(t) de- 
termines the maximum At by which the system may be 
propagated in a single time step. For time independent 
H , in principle, arbitrary large time steps are possible at 
the expense of increasing M. We may choose M such 
that for k > M the modulus of all expansion coefficients 
\c k (aAt/h)\ ~ J k (aAt/h) is smaller than a desired accu- 
racy cutoff. This is facilitated by the fast asymptotic decay 
of the Bessel functions, 



J k (aAt/h) 



for k — * oo 
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Thus, for large M, the Chebyshev expansion can be con- 
sidered as quasi-exact, and permits a considerably larger 
time step than e.g. the Crank-Nicholson scheme. Besides 
the high accuracy of the method, the linear scaling of com- 
putation time with both time step and Hilbert space di- 
mension are promising in view of potential applications to 
more complex systems. In our cases almost all computa- 
tion time is spent in sparse MVMs, which can be efficiently 
parallelised, allowing for a good speedup on parallel com- 
puters. 

2.3. Semiclassical methods 

During the last decades, a variety of semiclassical meth- 
ods have been tailored in order to incorporate certain quan- 
tum effects at least partially into classical many-particle 
simulations. Based on the real time (Feynman-) path in- 
tegral formulation of quantum mechanics, where action 
integrals take the center stage, they allow propagation of 
the (complex) wave function of a quantum system in time. 
Within the numerical evaluation of the integrals occur- 
ring by Monte Carlo (MC) techniques [13], the oscillatory 
complex valued integrand causes a dynamical sign problem 
which spoils the efficiency of the MC integration. 

Wigner-Moyal approach. Since a quantum system can 
be described equivalently in terms of real valued quantum 
phase space distribution functions e.g., the Wigner 



function, the dynamical sign problem may be alleviated [11 



[l2j].Fl Starting from the von Neumann equation we may de- 
rive an evolution equation for the Wigner function W(q, p, t) 
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where F(q) 
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is the classical force, and 
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(9) 

In the classical limit, the right hand side of (JSJ vanishes, 
leaving us with the Liouville equation for the phase space 
density. Then the dynamics can be expressed in terms of 
the classical propagator 



^ W (q,P,t;qo,Po,to) =S[q - q(t;p ,q ,t Q )} 

X S\p-p(t;p ,q Q ,t )] 



(10) 



where p and q are the momentum and coordinate of a tra- 
jectory that evolves according to the Hamilton equations 
of motion with initial conditions p(to) — Po and g(to) = qo ■ 



iNote that the Wigner function is just a convenient mathematical 
tool for the description of quantum systems and cannot be considered 
as a joint probability due to its possibly negative values, and conflict 
with Heisenberg uncertainty relation. 



Usin g II ^ , we may rewrite ((8]) in form of an integral equa- 
tion Bum, 

W(q,p,t) = J dpodqo U W (q,p,t]q ,p ,t )W {qo,Po,to) 
t 

+ J dr J dp T dq T U w (q,p,t;q T7 p T7 r) 
to 

00 

x / dsW(q T ,p T - s,T)ui(s,q T ) , 



(11) 

and solve it by iteration. Here, we consider only the lowest 
order, which means that the second integral is neglected 
completely. That is, we propagate classical trajectories 
(q,p) in time, after sampling their initial conditions po 
and qo from the initial Wigner function Wo(qo,Po, to) at 
time to by a MC procedure. Their superposition at the 
next time grid point gives W(q,p,t). The importance of 
higher order terms in the iteration series was investigated 
in Refs. [Hill- 

Quantum tomographic approach. As the dynamical 
sign problem is still present for the Wigner function, some 
years ago the description of quantum states in terms of 
a strictly positive function, the so called quantum tomo- 
gram, has been proposed [la]. Such a description seems 
promising in view of an effective MC sampling of the tra- 
[Fy]ectories during the propagation. The quantum tomo- 
gram @, m, 
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relates to the Wigner function by a class of Radon trans- 
forms [6] which arc characterised by \i and v. Each to- 
mogram contains a density information, and tuning (/z, v) 
appropriately, we may continuously switch between coordi- 
nate and momentum representation. Also for the quantum 
tomogram, an evolution equation can be derived from the 
von Neumann equation 
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For harmonic potentials, Eq. (| 13[) can be reformulated 
as a continuity equation and solved by trajectory meth- 
ods. For potentials of arbitrary shape, the complicated 
structure prevents a direct evaluation of (|13jl in order to 
get w(X, p., v, t). A possible way out is a local expansion 
of the potential up to second order [l9j]. The identifica- 
tion with the continuity equation then holds locally for 
each coordinate about which the potential is expanded. 
Since the slope and local curvature of the potential enter 
the propagator for the trajectories, the efficiency of the 
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tomogram-reconstruction depends crucially on the choice 
of the potential sampling positions. An intuitive, albeit 
not unique, choice for those positions are the coordinates of 
classically evolving particles similar to those in the Wigner 
approach. Advantageously, it is not necessary to propa- 
gate the whole set of tomograms if one is interested in the 
tomogram in one reference frame (/x, v) only. Instead, it is 
sufficient to find those trajectories which end up in a given 
(X (t) , n(t) , v(t)), or equivalently to propagate the desired 
(X, /i, ^-trajectories backward in time: one needs to be 
aware of the non-uniqueness of the propagator due to the 
various possible choices of the potential sampling points. 

3. Topical applications in physics 

We now apply the numerical techniques presented in 
the preceding section to selected physical systems and sit- 
uations. As a first step, we calibrate the different ap- 
proaches by studying a simple double well toy model. De- 
tecting the limitations and prospects of the various meth- 
ods seems to be necessary before applying them to the 
more complicated problems of current interest. Let us 
point out that the considered implementations of the semi- 
classical approaches provide an approximate description 
of quantum mechanics only, i.e. they will clearly not re- 
produce all the quantum effects. For the finite graphene 
nanoribbons studied in the second example it is barely pos- 
sible to determine the full quantum dynamics by means of 
direct diagonalisation techniques or applying the Crank- 
Nicholson scheme. This gives us the opportunity to bench- 
mark their performance in comparison to the Chebyshev 
expansion for a system for which the Hamilton matrix 
is not tridiagonal. The last example has been chosen to 
demonstrate the applicability of the Chebyshev approach 
to a true many-particle problem, the tunnelling of a po- 
laron. There the Hilbert space dimension is so large that 
neither the direct diagonalisation method nor the Crank- 
Nicholson can be applied. 

3.1. Double well potential 

As a basic test case we consider the motion of a Gaus- 
sian wave packet in the ID double well potential 

V(q) = V + l -mul{-q 2 + aq 4 ) , (14) 

sketched in the top panel of Fig. [1] We use m = u m , 
ui^Ut = 0.4, auj = 0.02 and Vq = ue = u m uj/u^, where 
u m , u t , U£ are the reference units for mass, time and length, 
respectively. The initial Gaussian of width a is centred at 
go, with centre-of-mass momentum po, 
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where we choose po — u m ui/ut, qo = — 5u£, and a = 
U(/\/2. The top panel of Fig. [T] gives also the real and 



Figure 1: Time evolution of a Gaussian wave packet in a double 
well potential (top panel). Results displayed are obtained using the 
Chebyshev, the first order Wigner-Moyal and the tomographic ap- 
proach (from top to bottom), respectively. 
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imaginary part of the wave function, where the baseline 
indicates the energy of the initial state in relation to V(q). 

Discretising the potential V(q) on an equally spaced 
grid of N = 2048 q-points, we fully diagonalise the Hamil- 
tonian to get the exact dynamics as reference. If we choose 
the iterative time step accordingly, the real and imaginary 
part of the wave function at t — AQu t is reproduced by both 
iterative methods with an absolute error less than 10 -9 . 
For this the maximum allowable time step for the Crank- 
Nicholson scheme is 5t = 4 x lCP 6 itt. For the Cheby- 
shev approach, the accuracy is even better than 2 x 10~ n . 
Here the required time step is related to the order of the 
Chebyshev expansion, i.e. the number of moments M. 
For M = 1500 we may propagate the wave function by 
At = 0.4w t . At the expense of calculating a larger number 
of moments M, larger times steps may be chosen without 
loss of accuracy. This will be demonstrated in the next 
section. The time evolution of the modulus squared of the 
wave function is shown in the second panel of Fig.[TJ While 
the major part of the wave packet stays in the left well, 
a sizeable fraction also penetrates the barrier. The rich 
structure in the density pattern reflects the (well-known) 
presence of strong interference effects. 

Restricting the Wigner-Moyal scheme to first order, the 
corresponding data reproduces in essence the overall dy- 
namics, but the fine interference patterns are not resolved 
within this approximation (see third panel of Fig. [T]). There 
are two major parameters which influence the computa- 
tion time for this approach: (a) the number of propagated 
trajectories and (b) the time step necessary for their clas- 
sical propagation. For the results presented, we have used 
5 x 10 5 trajectories and a propagation time step of 0.04it t . 
The agreement between the exact solution and the tomo- 
graphic result (see also the bottom panel of Fig. QJ is even 
worse as the tunnelling to the right hand side of the barrier 
is missed almost completely. Furthermore, the probability 
density for large negative values is overestimated, i.e. the 
effect of the steep anharmonic confinement potential is not 
accounted for correctly. In Tab. [T] we summarise the run 
times T run on a single Xeon 5160 processor required by the 
different methods in order to follow the time evolution of 
the system up to t = 40to- 

This very basic example already shows that a straight- 
forward use of the Wigner-Moyal and tomographic ap- 
proaches only partially accounts for quantum effects. While, 
in principle, both methods are equivalent with respect to 
the solution of the time dependent Schrodinger equation, 
an efficient implementation is lacking. For the Wigner- 
Moyal formalism there are two prospects. If the accuracy 
of the first order approximation is satisfactory, i.e. the 
neglect of the fine interference patterns is tolerable, this 
method provides an acceptable performance and might 
show its true virtue for many-particle systems. If one 
has to include higher terms of the iteration series, how- 
ever, e.g. because subtle quantum effects are important, 
the method is not competitive anymore because (i) the 
computational requirements increase drastically and (ii) 
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Table 1: Time evolution of a wave packet in the double well potential 
up to t = 40to- Data gives the run times by exact diagonalisation 
(ED), Crank-Nicholson (CN), Chebyshev (C), Wigner-Moyal (WM) 
and tomographic (T) methods. 



numerical fluctuations amplify strongly during the calcu- 
lation [19( | . The practical applicability of the tomographic 
approach to arbitrarily shaped potentials is also question- 
able, mainly because there is no simple way to construct 
suitable sampling functions for the coordinate sampling 
in the potential evaluation. Extensions beyond harmonic 
potentials will suffer crucially from this limitation. Apart 
from the poor accuracy of the results, also the compu- 
tational requirements were significant higher than for the 
other methods (although we used only 800 trajectories in 
this calculation). 

3.2. Disordered graphene nanoribbons 

Recently much interest has been devoted to investi- 
gate how disorder influences the transport properties of 
graphene [It], H3- ^ * s known that the presence of ar- 
bitrarily weak disorder leads to Anderson localisation of 
the single particle wave function on infinite 2D square lat- 
tices |l| . In weakly disordered 2D systems the localisation 
lengths are huge, however, and may easily become compa- 
rable to the system sizes that are technologically relevant 
e.g. for graphene nanoribbons. In those cases, we expect 
a conducting behaviour of the device despite the presence 
of disorder causing localisation on larger length scales. 

To investigate the influence of Anderson disorder for 
finite graphene nanoribbons, we consider a tight-binding 
model on a honeycomb lattice, 



H 



H.c. 



N 
3 = 1 



(16) 



Here the operators cj (cj) create (annihilate) an electron 
in a Wannier state centred at site j. The on-site poten- 
tials Cj are random variables in the interval [-W/2, W/2], 
where W is a measure for the disorder strength. The elec- 
tron transfer between nearest-neighbour lattice sites (ij) 
is described by the transfer integral t. Tailoring stripes 
(ribbons) out of the infinite honeycomb-lattice, we have to 
distinguish two cases with respect to the boundary condi- 
tions. Depending on the orientation of the stripe we get 
boundaries of either zigzag or armchair type. Since in ex- 
perimental probes armchair edges are more common, we 
will focus on those in the following. 

Starting from a wave packet which is localised on one 
site in the centre of each ribbon, we evolve the quantum 
state using the Chebyshev approach described in Sect. 12.21 
We consider devices of 1.11 x 212.8 nm 2 with 10 x 1000 
atoms. The main panel of Fig. [5] displays the time evo- 
lution of the wave packet for one realization of disorder 
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Figure 2: Time evolution of a localised state on an armchair graphcnc 
nanoribbon for different values of disorder. The left/right asymmetry 
of the state is due to the particular disorder realisation and initial 
state. 



for several values of disorder strength (time is measured in 
units of the inverse hopping element, to = 1/t)- 

The following aspects of the wave function dynamics 
should be noted: (i) The initially localised wave function 
spreads with time and reaches its maximum extension at 
about t ~ 10 3 io, independent of the disorder strength. 
Also for much longer times (t ~ 10 4 to) this extension 
does not change significantly anymore, (ii) The disorder 
strength strongly influences the spatial region over which 
the state spreads, i.e. the localisation length. While we see 
clear evidence for localisation at large disorder (W = 4i), 
the localisation length for weak disorder (W — 0.5?) is 
markedly larger than the system size, leading to an evenly 
spread state on the ribbon. Note that the occurrence of an 
extended (conducting) state in the latter case is only due 
to the finite system size. For longer ribbons a disorder of 
W = 0.5t is sufficient to localise the wave function as well. 

The disordered nanoribbon setup gives us a good op- 
portunity to benchmark the Chebyshev approach against 
the Crank-Nicholson scheme. Since the Hamiltonian is no 
longer tridiagonal the solution of the linear equation sys- 
tem cannot be done by the Thomas algorithm. Instead, we 
use the standard solver for double-complex linear equation 
systems from LAPACK, ZGESV. Table [U summarises the 
number of moments required to get agreement between the 
Chebyshev and the exact results for different time steps 
At. The computation times for calculating the quantum 
state at time t = 10 4 t using the various At is also given. 
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Table 2: Numeber of Chebyshev moments (M) and overall runtime 
(At /to) required for the calculation of the time evolution up to t = 
10 4 to on a single Intel Xeon 5160 core. 



For comparison, we also give run time T run for the exact 
diagonalisation. As one iteration of the Crank-Nicholson 
scheme using ZGESV takes 4 minutes, the given T mn is 
only an estimate. 

3.3. Spatiotemporal evolution of polaron states in finite 
quantum structures 

The cycle of quasiparticle formation is fundamental to 
many fields of physics. In condensed matter, e.g., the 
coupling between a charge carrier and the lattice degrees 
of freedom may create a new quasiparticle, an electron 
dressed by an phonon cloud. This composite entity is 
called polaron. From the basic electron-phonon (EP) in- 
teraction processes, the absorption/emission of a phonon 
with a simultaneous change of the electron state, it is clear 
that the motion of even a single electron in a deformable 
lattice constitutes a complex many-body problem, in that 
phonons are excited at various positions, with highly non- 
trivial dynamics [Til ]. Polaron transport through finite 
quantum systems becomes increasingly important for nan- 
otechnology applications. 

The microscopic structure of polarons is very rich. Fo- 
cusing on polaron formation in systems with short-range 
non-polar EP interaction and site-dependent potentials, 
we consider a generalised Holstein Hamiltonian [j| 

H = ^Aini-t^cJcj+i+Hx.) 

i i 

-^9iUa{b\ + b^nie + u ^2 bjbi , (17) 

i.o i 

where c\ (b\) creates an electron (phonon) at site i of a 
ID lattice. Parameters are the electron transfer integral 
i, the EP coupling strength j?j = [(e p + Ep^/wo] 1 / 2 , and 
the phonon frequency loq. The potentials Aj can describe 
a tunnel barrier, disorder, or a voltage basis. 

How does a bare electron time evolve to become a po- 
laron quasiparticle? To what extent can a polaron tunnel 
through a quantum barrier? Having the iterative Chebyshev- 
based time evolution algorithm explained in Sec. 2.2 at 
hand, we can address these questions in the framework 
of model (|17p . Let us emphasise that our numerical ap- 
proach, acting in the tensorial product Hilbert space of 
electron and phonons, takes into account the full dynam- 
ics of both quantum objects. Since the Hilbert space as- 
sociated to the phonons is infinite, we applied a controlled 
truncation procedure retaining only basis states with at 
most iVph phonons ID, . However the truncated Hilbert 
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space dimension (D tot ) is still very large even for small lat- 
tices and the dimension of the corresponding sparse matrix 
problem does limit the physical parameter region attain- 
able. Thus, we use a memory saving implementation of the 
sparse MVM where the non-zero matrix elements are not 
stored but recomputed in each sparse MVM step, limiting 
the overall memory consumption of our implementation to 
five vectors of size D tot . In this context we can access a 
massively parallel sparse MVM code which has proven to 
be sufficient to compute the ground state of the model (fTTf 
up to D tot = 3.5 x 10 11 very efficiently on more than 5000 
processor cores jjj. For the single polaron dynamics pre- 
sented here, the matrix dimension is D tot — 6.2 x 10 8 
and we run the Chebyshev approach on 18 processors of 
a SGI Altix4700 compute server accessing a total of ap- 
proximately 60 GBytes of main memory and consuming 
less than 500 CPU-hrs to compute the result presented in 
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Figure 3: Quantum dynamics of polaron formation and polaron tun- 
nelling through a potential barrier. A Gaussian wave packet cen- 
tred around site 4 begins evolving at time t = with momentum 
k = it/9. Moving to the right, a polaron is formed that hits the 
barrier located at site 12 at about t ~ 10. There complicated reflec- 
tion and transmission processes of the composite quasiparticle made 
up of an electron and phonons take place. Displayed is the time- 
evolution of the local particle densities (rij) (solid lines, open circles), 
phonon numbers (b|6 i ) (dot-dashed lines, stars) and EP correlations 
yff? = (n-(bi + bi~)) (dashed lines). Open boundary conditions were 
used at sites 1 and 18. In the numerics we account for all states with 
up to iVph < 11 phonons and in the ground state the weight of all 
basis states containing exactly N p ^ = 11 is less than 10 -11 . 

Figure [3] shows snapshots of polaron formation and po- 
laron propagation at intermediate EP coupling s p = 1 in 
the non-adiabatic regime luq = 2 (the time and all ener- 
gies were measured in units of f -1 and i, respectively). At 
t = a bare electron wave packet is injected at site 4 and 
launched to the right. Shortly after, the electron is not yet 
dressed and moves nearly as fast as a free particle. But 
then the electron emits (creates) phonons in the vicinity 



of the electron's starting point, in order to reduce its en- 
ergy to near the bottom of the band, and then forms a 
polaron (see the panel at t = 6). One of the most impor- 
tant properties of the polaron is an increased inertial mass, 
for the reason that some phonons have to travel with the 
particle (as indicated by the enhanced on-site EP corre- 
lations). At the same time we observe a "backscattered" 
current [l4T | , evolving to a left moving polaron. When the 
right-moving polaron reaches the wall at site 12 it will be 
partly reflected. More importantly, the additional local 
EP interaction £ p ,i2 renormalises the on-site adiabatic po- 
tential at site 12, i.e. leads to a local polaron level shift 
that softens the barrier A12 — 2. As a result a vibration- 
mediated tunnelling of the particle takes place, whereby 
some phonons are stripped at the barrier and are recol- 
lected by the transferred particle afterwards (cf. the snap- 
shots from t = 10 to 14). Finally, the particle is reflected 
at the boundary and moves to the left passing the bar- 
rier again. Note that during the whole run time uncorre- 
cted phonon excitations remain in the system, especially 
near the injection point. In our opinion this example im- 
pressively demonstrates that our approach can be used to 
monitor the complicated multiple time-scale dynamics of 
quasiparticle transport in finite quantum structures. 

4. Conclusion 

To summarise, in this work we compared various nu- 
merical approaches to the dynamics of complex quantum 
systems: expansion into eigenstates, iterative Crank-Nicholson 
and Chebyshev schemes, as well as semiclassical Wigner- 
Moyal and quantum tomography methods. The different 
methods have been applied to several physical systems and 
problems, ranging from motion in a simple double- well 
toy model to questions of current interest such as electron 
transport in disordered graphene nanoribbons or polaron 
motion in a finite quantum structure. 

The Wigner-Moyal approach, evaluated in first order 
of the iteration series, essentially reproduces the quan- 
tum dynamics. Nevertheless important quantum interfer- 
ence effects, appearing in the exact solution, are missed. 
The successful application of the quantum tomogram to 
the time evolution of quantum systems crucially depends 
on a suitable sampling algorithm for the coordinates at 
which the potential is evaluated. If those are sampled ac- 
cording to trajectories of classically propagated particles, 
the result for the quantum dynamics is rather poor and 
fails to describe tunnelling and anharmonicity effects cor- 
rectly. While the moderate numerical costs of the first or- 
der Wigner-Moyal approach seem appealing for a possible 
application to more complex systems, the computational 
resources required by the quantum tomographic approach 
are high in general. 

On the side of the exact techniques, the Chebyshev ap- 
proach largely outperforms the standard Crank- Nicholson 
scheme, both in computation speed (usable time step) and 
treatable system sizes (only matrix- vector multiplications 
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were required). We conclude that the Chebyshev approach 
represents a very efficient and reliable tool to determine 
the quantum dynamics even for rather complex interact- 
ing many-particle systems. 



References 

[1] E. Abrahams, P. W. Anderson, D. C. Licciardello, and 

T. V. Ramakrishnan. Phys. Rev. Lett, 42:673, 1979. 
[2] A. S. Arkhipov and Y. E. Lozovik. Phys. Lett. A, 319:217, 

2003. 

[3] A. S. Arkhipov, Y. E. Lozovik, and V. I. Man'ko. Phys. 

Lett. A, 328:419, 2004. 
[4] R. Chen and H. Guo. Comp. Phys. Comm., 119:19-31, 

1999. 

[5] J. K. Cullum and R. A. Willoughby. Lanczos Algo- 
rithms for Large Symmetric Eigenvalue Computations, 

Birkhauser, Boston, 1985. 

[6] S. R. Deans. The Radon transform and some of its appli- 
cations. John Wiley & Sons, New York, 1983. 

[7] H. Fehske, A. Alvermann, and G. Wellein. In S. Wagner, 

M. Steinmetz, A. Bode, and M. Brehm, editors, High Per- 
formance Computing in Science and Engineering, Garch- 
mg/Munich 2007, pages 649-668. Springer- Verlag, Berlin, 
2009. 

[8] H. Fehske, A. Weifie, and R. Schneider, editors. Computa- 
tional many-particle physics, volume 739 of Lecture Notes 

in Physics. Springer, Berlin Heidelberg, 2008. 
[9] H. Fehske, G. Wellem, J. Loos, and A. R. Bishop. Phys. 

Rev. B, 77:085117, 2008. 
10] V. S. Filinov. Nucl. Phys. B, 271:717, 1986. 
11 V. S. Filinov. Mol. Phys., 88:1517 & 1529, 1996. 
12] V. S. Filinov, Y. V. Medvedev, and V. L. Kamskyi. Mol. 

Phys., 85:711, 1995. 
[13] E. 'Jeckelmann and H. Fehske. Rivista del Nuovo Cimento, 

30:259, 2007. 

[14] L.-C. Ku and S. A. Trugman. Phys. Rev. B, 75:014307, 
2007. 

[15] H. W. Lee. Physics Reports, 259:147, 1995. 

[16] S. Mancini, V.l. Man'ko, and P. Tombesi. Phys. Lett. A, 

213:1, 1996 & Found. Phys., 27:801, 1997. 
[17] V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, 

N. M. R. Peres, and A. H. Castro Neto. Phys. Rev. Lett., 

96:036801, 2006. 
[18] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. 

Vetterling. Numerical recipes. Cambridge University Press, 

Cambridge, 1986. 
[19] G. Schubert, V. S. Filinov, K. Matyash, R. Schneider, and 

H. Fehske. arXiv:0810.4302 , preprint 2008. 
[20] H. Tal-Ezer and R. Kosloff. f. Chem. Phys., 81:3967-3971, 

1984. 

[21] A. Weifie and H. Fehske. Lecture Notes in Physics, 
739:545, 2008. 

[22] G. Wellein, H. Roder, and H. Fehske. Phys. Rev. B, 
53:9666, 1996. 

[23] S.-J. Xiong and Y. Xiong. Phys. Rev. B, 76:214204, 2007. 



8 



